function accuracy = cca_window_harmonics_select(Fs,stimulusFreq,nHarmonics,epochLengths,stepTime,data)
%画出不同窗长与谐波次数下的CCA正确率图，用来选择合适的窗长与谐波次数
%---input---
%Fs：采样率,单位Hz，e.g.500
%stimulusFreq：刺激频率，单位Hz，e.g.[7.5 9.0 10.5 12.0]
%nHarmonics:谐波次数，e.g. 1:8，目前不能超出8
%epochLengths：窗长范围，单位s，e.g.0.5:0.1:2.0
%stepTime：滑窗长度，单位s，e.g.0.2
%dataa：输入数据，第一个维度为session，第二个维度为stimulus，第三个维度为trial，第四个维度为channel，第五个维度为samplepoints
%---output---
%accuracy：正确率，size1窗长，size2谐波次数，size3刺激频率


for L = 1:length(epochLengths)
    epochTime = epochLengths(L);
    for h = 1:length(nHarmonics)
        result = [];
        for stimulus = 1:size(data,2)
            for session = 1:size(data,1)
                for trial = 1:size(data,3)
                    epoch = data(session,stimulus,trial,:,:);  %单一trial数据
                    
%                     %移除错误数据,根据数据的不同也不同
%                     if((session == 2) && (trial == 10) && (stimulus == 2) && (condition == 2))
%                         epoch = data(session,stimulus,trial-1,:,:);
%                     end
                    
                    epoch = reshape(epoch,size(data,4),size(data,5));
                    windowCount = floor((size(epoch,2) - epochTime * Fs + 1) / (stepTime * Fs));
                    parfor window = 1:windowCount
                        result(session,stimulus,trial,window) = canonical_correlation_analysis(epoch(:,int32((window-1) * stepTime * Fs + 1):int32((window-1) * stepTime * Fs + epochTime * Fs)),stimulusFreq,h,epochTime,Fs);
%                         display(['epochLength: ', num2str(epochTime), ' nHarmonics: ', num2str(h), ' Stimulus: ', num2str(stimulus), ...
%                             ' Session: ', num2str(session), ' Trial: ', num2str(trial), ' Window: ', num2str(window), ...
%                             ' Result: ', num2str(result(session,stimulus,trial,window))]);
                    end
                end
            end
        end
        %计算正确率
        for stimulus = 1:size(data,2)
            accuracy(L,h,stimulus) = sum(result(:,stimulus,:,:) == stimulus, 'all') / (size(result,1) * size(result,3) * size(result,4)) * 100;
        end
        display(['epochLength: ', num2str(epochTime), ' nHarmonics: ', num2str(h),' accuracy: ', num2str(mean(accuracy(L,h,:)))]);
    end
end

%画图
Marker = {'o', '+', '*', 'x', '^', '<', '>', 'v'};
legendString = cell(length(nHarmonics),1);
for n = 1:length(legendString)
    legendString{n} = strcat('Nh = ',num2str(nHarmonics(n)));
end
acc = mean(accuracy, 3);
acc = reshape(acc, length(epochLengths), []);
figure;
for nh = 1:length(nHarmonics)
    plot(epochLengths,acc(:,nh),'linewidth',2,'Marker',Marker{nh});
    hold on;
end
xlabel('Epoch Length [s]');
ylabel('Accuracy %');
legend(legendString);
set(gca, 'linewidth', 2, 'fontsize', 20, 'fontname', 'times');

end









